!// Dmitry.Konovalov@jcu.edu.au Dmitry.A.Konovalov@gmail.com, 2013

!
! This is for any general atom data/calcs
!
module jm_atom_mod;
use vec_mod;
implicit none;
logical, save :: ATOM_TRACE = .TRUE.;
private;

! NOTE 12-Apr-2013: name any new public functions as atom_xxx()
public ATOM_ELEC_SPIN2, atom_wfv_h1s, atom_wfv_hnL, atom_hlike_eng &
&, atom_z1v_1s, atom_z0v_1s, atom_y0v_1s, atom_y0v_2s, atom_y0v_2p, atom_y2v_2p;

integer, parameter :: ATOM_ELEC_SPIN2 = 1; ! spin*2 for electron



contains;!=====================================================================


! n - principal quantum number
! L - orbital angular momentum 
! nucZ - nuclear electric charge (eg '+1' for atomic hydrogen)
! orbZ - orbiting particle charge (eg '-1' for electron)
real(DP) function atom_hlike_eng(n, L, nucZ, orbZ) result(res);         
  real(DP), intent(IN) :: nucZ, orbZ;      
  integer, intent(IN) :: n, L;      
  real(DP) :: z, dn;
  call dbg_assert_true("L < n", L < n); !// 1s, 2s, 2p, 3s, 3p, 3d
  call dbg_assert_true("n > 0", n > 0);
  call dbg_assert_true("L .GE. 0", L .GE. 0);
  z = nucZ * orbZ;
  call dbg_assert_true("nucZ * orbZ < D0", z < D0);
  dn = real(n, DP);
  res = - D05 * z * z / (dn * dn); 
  return;
end function


SUBROUTINE atom_wfv_hnL(res, rv, n, L, absZ);     
  real(DP), intent(IN)  :: rv(:), absZ;      
  integer, intent(IN)  :: n, L;      
  real(DP), intent(OUT) :: res(:);  
  call dbg_assert_true("absZ > D0", absZ > D0);
  call dbg_assert_not_zero("absZ .ne. zero", absZ);
  if (n==1 .AND. L==0) then
    call atom_wfv_h1s(res, rv, absZ);
  elseif (n==2 .AND. L==0) then
    call atom_wfv_h2s(res, rv, absZ);
  elseif (n==2 .AND. L==1) then
    call atom_wfv_h2p(res, rv, absZ);
  else
    call dbg_int("n=", n);
    call dbg_int("L=", L);
    call dbg_assert_true("atom_wfv_hnL not ready", .FALSE.);
  endif
end SUBROUTINE



! wfv-wave function vector; h- hydrogen-like; 1s-bound state
! rv-radial vector
! absZ - abs-value if the hydrogen like nuclear charge 
! norm = \int_0^\infty wf^2 dr = 1; 
SUBROUTINE atom_wfv_h1s(res, rv, absZ);     
  real(DP), intent(IN)  :: rv(:), absZ;      
  real(DP), intent(OUT) :: res(:);  
  integer :: i, mini, maxi;
!  call dbg_str("-->atom_wfv_h1s()", ATOM_TRACE);
  call dbg_assert_true("absZ > D0", absZ > D0);
  call dbg_assert_not_zero("absZ .ne. zero", absZ);
  call vec_init_res(res, rv);
  call vec_range(mini, maxi, rv); 
  do i = mini, maxi 
    res(i) = atom_wf_h1s(rv(i), absZ);
  end do
end SUBROUTINE
! wf-wave function; h- hydrogen-like;
! norm = \int_0^\infty wf^2 dr = 1; 
real(DP) function atom_wf_h1s(r, absZ) result(y);         
  real(DP), intent(IN) :: r, absZ;      
  real(DP) :: rz;      
  rz = r * absZ;
  y = D2 * SQRT(absZ) * EXP(-rz) * rz;
  return;
end function

SUBROUTINE atom_wfv_h2s(res, rv, absZ);     
  real(DP), intent(IN)  :: rv(:), absZ;      
  real(DP), intent(OUT) :: res(:);  
  integer :: i, mini, maxi;
!  call dbg_str("-->atom_wfv_h1s()", ATOM_TRACE);
  call dbg_assert_true("absZ > D0", absZ > D0);
  call dbg_assert_not_zero("absZ .ne. zero", absZ);
  call vec_init_res(res, rv);
  call vec_range(mini, maxi, rv); 
  do i = mini, maxi 
    res(i) = atom_wf_h2s(rv(i), absZ);
  end do
end SUBROUTINE
! wf-wave function; h- hydrogen-like;
! norm = \int_0^\infty wf^2 dr = 1; 
real(DP) function atom_wf_h2s(r, absZ) result(y);         
  real(DP), intent(IN) :: r, absZ;      
  real(DP) :: rz;      
  rz = r * absZ * D05;
!  double d = r * z * 0.5;
!  return Math.sqrt(2 * z) * Math.exp(-d) * d * (1.0 - d);
  y = SQRT(D2 * absZ) * EXP(-rz) * rz * (D1 - rz);
  return;
end function



SUBROUTINE atom_wfv_h2p(res, rv, absZ);     
  real(DP), intent(IN)  :: rv(:), absZ;      
  real(DP), intent(OUT) :: res(:);  
  integer :: i, mini, maxi;
!  call dbg_str("-->atom_wfv_h1s()", ATOM_TRACE);
  call dbg_assert_true("absZ > D0", absZ > D0);
  call dbg_assert_not_zero("absZ .ne. zero", absZ);
  call vec_init_res(res, rv);
  call vec_range(mini, maxi, rv); 
  do i = mini, maxi 
    res(i) = atom_wf_h2p(rv(i), absZ);
  end do
end SUBROUTINE
real(DP) function atom_wf_h2p(r, absZ) result(y);         
  real(DP), intent(IN) :: r, absZ;      
  real(DP) :: rz;      
  rz = r * absZ * D05;
!  double d = r * z * 0.5;
!  return Math.sqrt(2 * z / 3.) * Math.exp(-d) * d * d;
  y = SQRT(D2 * absZ / D3) * EXP(-rz) * rz * rz;
  return;
end function


SUBROUTINE atom_y0v_1s(res, rv);     
  real(DP), intent(IN)  :: rv(:);      
  real(DP), intent(OUT) :: res(:);  
  integer :: i, mini, maxi;
  !call dbg_str("-->atom_wfv_h1s()", ATOM_TRACE);
  call vec_init_res(res, rv);
  call vec_range(mini, maxi, rv); 
  do i = mini, maxi 
    res(i) = atom_y0_1s(rv(i));
  end do
end SUBROUTINE
!// Y^0(1s, 1s; r) = 1-exp(-2r)*(1+r)
real(DP) function atom_y0_1s(r) result(y);         
  real(DP), intent(IN) :: r;      
  y = D0;
  if (dbg_is_zero(r) .OR. r < D0) return;     
!  return 1. - Math.exp(-2. * r) * (1. + r);
  y = D1 - EXP(-D2 * r) * (D1 + r);
  return;
end function


SUBROUTINE atom_y0v_2s(res, rv);     
  real(DP), intent(IN)  :: rv(:);      
  real(DP), intent(OUT) :: res(:);  
  integer :: i, mini, maxi;
  !call dbg_str("-->atom_wfv_h1s()", ATOM_TRACE);
  call vec_init_res(res, rv);
  call vec_range(mini, maxi, rv); 
  do i = mini, maxi 
    res(i) = atom_y0_2s(rv(i));
  end do
end SUBROUTINE
!// Y^0(1s, 1s; r) = 1-exp(-2r)*(1+r)
real(DP) function atom_y0_2s(r) result(y);         
  real(DP), intent(IN) :: r; 
  real(DP) :: r3;     
  y = D0;
  if (dbg_is_zero(r) .OR. r < D0) return;     
  r3 = r**3;
!  //return 1. - exp(-r) * (1. + 3./4 * r + 1./4 * r * r + 1./8 * pow(r, 3));
!  return 1. - Math.exp(-r) * (1. + 3. / 4 * r + 1. / 4 * r * r + 1. / 8 * Mathx.pow(r, 3));
  y = D1 - EXP(-r) * (D1 + D3 / D4 * r + D1 / D4 * r * r + D1 / D8 * r3);
  return;
end function


SUBROUTINE atom_y0v_2p(res, rv);     
  real(DP), intent(IN)  :: rv(:);      
  real(DP), intent(OUT) :: res(:);  
  integer :: i, mini, maxi;
  call vec_init_res(res, rv);
  call vec_range(mini, maxi, rv); 
  do i = mini, maxi 
    res(i) = atom_y0_2p(rv(i));
  end do
end SUBROUTINE
real(DP) function atom_y0_2p(r) result(y);         
  real(DP), intent(IN) :: r; 
  real(DP) :: r3;     
  y = D0;
  if (dbg_is_zero(r) .OR. r < D0) return;    
!  double er = Math.exp(r);
!  // exp(x) = 1 + x/1! + x^2/2!
!  // (exp(x)-1)/x
!  //  = 1/1! + x/2! + x^2/3! + x^3/4! + x^4/5!  + x^5/6!  + x^6/7! + x^7/8!
!  //  = 1    + x/2  + x^2/6  + x^3/24 + x^4/120 + x^5/720 + x^6/5040

 
  r3 = r**3;
!  //return 1. - exp(-r) * (1. + 3./4 * r + 1./4 * r * r + 1./24 * pow(r, 3));
!  return 1. - Math.exp(-r) * (1. + 3. / 4 * r + 1. / 4 * r * r + 1. / 24 * Mathx.pow(r, 3));
  y = D1 - EXP(-r) * (D1 + D3 / D4 * r + D1 / D4 * r * r + D1 / D24 * r3);
  
  return;
end function


SUBROUTINE atom_y2v_2p(res, rv);     
  real(DP), intent(OUT) :: res(:);  
  real(DP), intent(IN)  :: rv(:);      
  integer :: i, mini, maxi;
  character(LEN=64) :: name = "atom_y2v_2p";
  call dbg_in(name);
  call vec_init_res(res, rv);
  call vec_range(mini, maxi, rv); 
  do i = mini, maxi 
    res(i) = atom_y2_2p(rv(i));
  end do
end SUBROUTINE
real(DP) function atom_y2_2p(r) result(y);         
  real(DP), intent(IN) :: r; 
  real(DP) :: r3, er, ex, c4, c5, c6, c7, c8, c9, c10, c11, c12, a, b;     
  y = D0;
  !call dbg_dp("atom_y2_2p(r=", r);
  if (dbg_is_zero(r) .OR. r < D0) return;     
  if (r < D01) then
    er = EXP(r);    call dbg_dp("er = EXP(r)=", er);
    !  // exp(x) = 1 + x/1! + x^2/2!
    !  // (exp(x)-1)/x
    !  //  = 1/1! + x/2! + x^2/3! + x^3/4! + x^4/5!  + x^5/6!  + x^6/7! + x^7/8!
    !  //  = 1    + x/2  + x^2/6  + x^3/24 + x^4/120 + x^5/720 + x^6/5040
    r3 = r * r * r;         !call dbg_dp("r3=", r3);
    c4 = D1/D120 - D1/D144; !call dbg_dp("c4=", c4);
    c5 = r / D720;
    c6 = c5 * r / D7;
    c7 = c6 * r / D8;
    c8 = c7 * r / D9;
    c9 = c8 * r / D10;
    c10 = c9 * r / D11;
    c11 = c10 * r / D12;
    c12 = c11 * r / D13;  !call dbg_dp("c12=", c12);
    b = r3 * (c4 + c5 + c6 + c7 + c8 + c9 + c10 + c11 + c12) / er; !call dbg_dp("b=", b);
    y = D30 * b;
  else 
    er = EXP(-r);    call dbg_dp("er = EXP(-r)=", er);
    !// exp(x) = 1 + x/1! + x^2/2!
    ex = (D1 - er)/r; !//
    !// (exp(x)-1)/x
    !//  = 1/1! + x/2! + x^2/3! + x^3/4! + x^4/5!  + x^5/6!  + x^6/7!
    !//  = 1    + x/2  + x^2/6  + x^3/24 + x^4/120 + x^5/720 + x^6/5040

    a = (D1 + (D1 / D2 + (D1 / D6 + (D1 / D24  + D1 / D144 * r) * r) * r) * r);
    b = ex - a * er;
    y = D30* (b / r);

    !//  double ex = (Math.exp(-r) - 1.)/(-r); // (exp(-r) - 1)/(-r)
    !//double a = (1.
    !//+ (1. / 2.
    !//+ (1. / 6.
    !//+ (1. / 24.  + 1. / 144. * r  ) * r
    !//) * r
    !//) * r
    !//);
    !//double b = ex - Math.exp(-r) * a;
    !//double res = 30.* (b / r);
  endif
  return;
end function


SUBROUTINE atom_z0v_1s(res, rv);     
  real(DP), intent(IN)  :: rv(:);      
  real(DP), intent(OUT) :: res(:);  
  integer :: i, mini, maxi;
  !call dbg_str("-->atom_wfv_h1s()", ATOM_TRACE);
  call vec_init_res(res, rv);
  call vec_range(mini, maxi, rv); 
  do i = mini, maxi 
    res(i) = atom_z0_1s(rv(i));
  end do
end SUBROUTINE
!// = INTL_0^r 4 s^2 exp(-2s) ds
!// INTL r^2 exp(ax) dx = exp(ax) (r^2/a - 2x/a^2 + 2/a^3)
!// a = -2
real(DP) function atom_z0_1s(r) result(y);         
  real(DP), intent(IN) :: r;      
  real(DP) :: a, ex; 
  y = D0;
  if (dbg_is_zero(r) .OR. r < D0) return;     
  a = -D2;
  ex = EXP(a * r);
  !return 4. * ex * (r * r / a - 2. * r / (a * a)) + (ex - 1) * 8. / (a * a * a);
  y = D4 * ex * (r * r / a - D2 * r / (a * a)) + (ex - D1) * D8 / (a * a * a);
end function


SUBROUTINE atom_z1v_1s(res, rv);     
  real(DP), intent(IN)  :: rv(:);      
  real(DP), intent(OUT) :: res(:);  
  integer :: i, mini, maxi;
!  call dbg_str("-->atom_wfv_h1s()", ATOM_TRACE);
  call vec_init_res(res, rv);
  call vec_range(mini, maxi, rv); 
  do i = mini, maxi 
    res(i) = atom_z1_1s(rv(i));
  end do
end SUBROUTINE
!// = INTL_0^r 4 s/r * s^2 exp(-2s) ds
!// INTL x^3 exp(ax) dx = exp(ax)/a (x^3 - 3x^2/a + 6x/a^2 - 6/a^3)
!// a = -2
real(DP) function atom_z1_1s(r) result(y);         
  real(DP), intent(IN) :: r;      
  real(DP) :: a, ex, ra4; 
  y = D0;
  if (dbg_is_zero(r) .OR. r < D0) return;     
  a = -D2;
  ex = EXP(a * r);
  ra4 = D24 / (r * a**4);
  y = D4 * ex / a * (r * r - D3 * r / a + D6 / (a * a));
  y = y + (D1 - ex) * ra4;
end function

end module